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Abstract 

The solution of the Lippman-Schwinger (L-S) integral equation is equivalent to the the solution 
of the Schrodinger equation. A new numerical algorithm for solving the L-S equation is described 
in simple terms, and its high accuracy is confirmed for several physical situations. They are: the 
scattering of an electron from a static hydrogen atom in the presence of exchange, the scattering 
of two atoms at ultra low temperatures, and barrier penetration in the presence of a resonance for 
a Morse potential. A key ingredient of the method is to divide the radial range into partitions, 
and in each partition expand the solution of the L-S equation into a set of Chebyshev polynomials. 
The expansion is called "spectral" because it converges rapidly to high accuracy. Properties of the 
Chebyshev expansion, such as rapid convergence, are illustrated by means of a simple example. 



1 



I. INTRODUCTION 



As stated in the textbook by Cummings, Laws, Redish and Cooney [jj, "Physics is a 
process of learning about the physical world by finding ways to make sense of what we 

n 

observe and measure. As the inspiring teacher Richard Feynman wrote, \2\ "Progress in all 
of the natural sciences depends on this interaction between experiment and theory"." 

An important tool required for carrying out this interaction is the solution of equations 
provided by a particular theory, in order to be able to compare its predictions with experi- 
ment. As the equations become more and more involved, such as in global climate study, in 
the construction of pharmaceutical drugs, in the analysis of large organic chains that exist 
in live cells, in the understanding of superconductivity, in the tracing of the earth's interior 
by means of seismic waves, in the construction of devices that transmit digital information, 
in the study of atomic, nuclear and particle theory (particularly in lattice gauge theory), 
etc., the resort to numerical computational methods becomes increasingly more necessary. 

The purpose of this paper is to point out special physical situations that require very 
accurate numerical algorithms, and to describe one such algorithm that has been recently 
developed. These special cases require either the evaluation of the solution of a wave equation 
out to large distances, or require high accuracy even for small distances, or both. Examples 
are the collision between atoms at extremely low temperatures. The understanding of such 
collisions is important for astro-physical applications, for the description of the state of atoms 
or molecules called Bose-Einstein condensates, and for the understanding of superfluidity in 
liquids formed out of weakly interacting atoms, such as the atoms of Helium. Helium is a 
"noble gas", i.e., its atoms interact mainly repulsively at short distances, yet, at intermediate 
distances (between 5 and 200 atomic units of distance) there is a small attractive valley in 
the potential energy curve (of a depth less than 3.5 x 10^^ atomic units of energy) within 
which a bound state can form. That weak attraction is in turn important for the molecular 
binding of a system of three or more helium atoms . The quantum mechanical wave 

function for the di-atom, in view of the weak binding energy of 4.4 x 10~^ atomic units of 
energy , extends to such large distances that accurate numerical values out to 2000 atomic 
units are required. 

For the case of the radial, one-dimensional, Schrodinger equation 

{(f/dr'^ + k'^)^l) = V%l), (1) 
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where k is the wave number in units of inverse length and V{r) is the potential in units 
of inverse length squared which contains the L{L + l)/r^ singularity, the most suitable 
equivalent integral equation for the S-IEM method is the Lippman-Schwinger equation 

^jJ{r) = sm{kr) + I go{r, r') V{r') tp{r') dr', (2) 







where Qq is the undistorted Green's function. In configuration space Qq has the well known 
semi-separable form Qo = —(l/k) sin(A;r<) cos(A;r>). (for negative energies one would have 
— (I/k) sinh(Kr<) exp(— Kr>)). By introducing the integral operator )Ct, so that when ap- 
plied on a function ip{r) the result is 

1 r I 

}CtiP{t) = —— cos{kr) I dr' sin(/cr') \^(r')V^(r') — — sin(A;r) / dr' cos(/cr') V{r')il){r'), 
k Jo k Jr 

(3) 

then Eq. Q can be written as 

■ip{r) = sin(A;r) + JCt "ipir), (4) 

where Kt ip means that ip{r') is included in the integrands contained in Eq. Q. This form of 
Eq. Q leads to the boundary condition that = 0, and since it assumes that for r >T 
the potential V{r) = 0, it leads to the asymptotic behavior ip{r) = sin{kr) + B cos{kr), 
where B is a constant determined from the solution of Eq. 0. If V{r) 7^ for r > T, 
then matching at r = T to the corresponding longrange functions (Bessel or Coulomb, for 
example) is required, as is explained in Ref. (Q], p])- 

A new method for solving the Lippman-Schwinger integral equation_P)). associated with 
the differential Schrodinger equation (^, has been developed recently jsl as an extension of 
a method due to Greengard and Rokhlin This method, to be called lEM (for integral 
equation method) has an accuracy which, for the same number of mesh-points, is far superior 
to the accuracy provided by finite difference methods for solving either an integral or a 

n 

differential equation. One of the intended applications [9] is the solution of the Faddeev 
equations for a three-body system in configuration space, since it requires the calculation 
of wave functions out to large distances. It is the purpose of this paper to describe the 
application of this method for positive energy, two-body scattering cases, and compare it 
with several other methods. The application of this method to finding bound-state negative 
energies is being developed, with the intention of obtaining the He-He bound state described 
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above. The basic idea of the lEM is to divide the radial interval into partitions, obtain two 
special solutions of the restricted Lippman-Schwinger equation in each partition, called 
Y{r) and Z{r), by expanding these solutions into a set of Chebyshev polynomials, and 
calculating the coefficients of the expansion in each partition. That expansion is " spectral" , 
i.e., it converges rapidly once the number of terms exceeds a certain value, and the error of 
truncating the expansion beyond that value is known, as is further explained below. Once 
the functions Y and Z are obtained in each partition, then the global function ip in that 
partition is expressed as a linear combination of the Y and Z . The coefficients of that 
combination are subsequently calculated by solving a matrix equation, which is sparse, as 
will be explained. Spectral expansions to solve integral equations, albeit using a rather 
different set-up, in particular not using Green's functions or partitions, has also recently 
been developed by B. Mihaila [l^ . 

Even though it is known that the errors which arise in the numerical solution of an 
integral equation are smaller than the errors in the solution of an equivalent differential 
equation, it is customary to solve the latter. The reason is that the algorithms for solving 
a differential equation by means of finite difference methods (such as Numerov of Runge- 
Kutta) are simple and do not require extensive storage space. By contrast, the discretization 
of an integral equation usually leads to large non-sparse matrices, and hence requires large 
investments of computer time and storage space. Therefore the gain in accuracy of the 
integral equation formulation is normally offset by a manifold increase in computational 
time. Our method circumvents this problem, as is described below. Before applications to 
physical cases are described, it is instructive to understand the basic accuracy properties of 
the spectral expansion method, as well as the basic ingredients of the lEM. 



II. SPECTRAL EXPANSION 



The main feature of a spectral expansion, namely its rapid convergence, will now be 
demonstrated by means of a simple example even though extensive discussions exist in the 



literature For the spectral expansion functions we will use Chebyshev polynomials 

only, although other orthogonal polynomials, such as Legendre, are also often used. We 
use Chebyshev polynomials because they are particularly well suited for obtaining the an- 
tiderivaties that appear in Eq. Q. 
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Spectral accuracy is described as follows: If a function f{x), — 1 < x < 1 is expanded 
in terms of Chebyshev polynomials Tj{x), 

oo 

m = ^ + J2a,T,{x) (5) 

then the error in truncating the expansion after n terms is proportional to + where p 
is the number of continuous derivatives which the function / has in the in the interval — 1 < 
X < 1. Furthermore, this truncation error is also proportional to the {n + l)'th coefficient 
of the expansion, which means that, after a certain number of terms, the coefficients aj 
decrease rapidly with j according to the same law In particular, if f{x) is infinitely 
differentiable, then the coefficients converge to zero asymptotically faster than any fixed 
power of Hence the term "spectral convergence" is also referred to as "superalgebraic 

convergence" . 

These properties will now be illustrated by expanding the function f{x) = exp(a;) into 
Chebyshev polynomials. The coefficients aj in Eq. ^ are given by 

2 f ^„/N/_ 0n_1/9, 2 



TT 



/ e"r,(x)(l-x')-^/'rfx= - r e^°^Vos(j^)rf^, (6) 

J-l TT Jo 



which follows from the orthogonality relation 



j Tk{x) Tj{x) (1 - x^Y^^ dx = ifj^k 



n/2 if j = k^0 

TT if j = k = 0. (7) 
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The integral in Eq. © can be calculated analytically. In view of Eq. (6.9.19) in Ref 
the result is aj = 2/^(1), where Ij{z) is a modified Bessel function of order j. Using the 
asymptotic expansion for large orders of a Bessel function, Eq. (9.3.1) of Ref. [lj|, an 
approximation to aj for large values of the index j is 

«i - 5 J ^ oo (8) 



Equation (jH)) shows that the value of aj decreases with j faster than any fixed power of j, as 
is also demonstrated in the Table HI . The first row lists the values of aj for j = 2, 4, 6, 8 as 
calculated from Eq. ©, (the results for the odd values of j are not shown) and the second 
row gives the values obtained from the asymptotic approximation (|H)) The table shows that 





«2 


a4 


ae 


as 




2.715S - 1 
2.60^ - 1 


5.474^ - 3 
5.32^ - 3 


4A50E - 5 
4.40^; - 5 


1.992^ - 7 
1.958^ - 7 



TABLE I: Chebyshev expansion coefficients a(k) of f(x)=exp(x) 
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FIG. 1: Truncation errors in the expansion of f{x) = exp{x) into Chebyshev polynomials, divided 
by the first expansion coefficient not included in the sum. 

the coefficients decrease rapidly with the order j. Will the truncation error also decrease 
rapidly? 

The truncation error in the expansion is defined as e„(a;) = f{x) — fn{x) where fn{x) 
denotes the sum in Eq. (0) that is taken from j = 1 to jmax = n — 1. A useful property 
of spectral expansions is that this error decreases with n proportionally to a„, the first 
expansion coefficient not included in the sum. This is demonstrated in Fig. ^ which shows 
the ratio en(x)/a„, for n = 2, 4, 6 and 8. The figure shows that the curves are approximately 
contained between ±1, i.e., the truncation error is of the same magnitude as a„ independently 
of the value of x. Hence the truncation error does not show a Gibbs phenomenon at the end 
points, as would be the case for an expansion into a Fourier Series. 

The above mentioned relation between the truncation error and the value of the Cheby- 
shev coefficient provides a convenient method for finding the appropriate size of each parti- 
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tion, compatible with the overall prescribed error. Clenshaw and Curtis who originated 
this spectral integration technique, recommend using the average size of the three last con- 
secutive coefficients as an accuracy criterion. 

Once the coefficients of the expansion (0) are known for j = 0, 1, ..A^, then one has a 
semi-analytical approximation to the function f{x), given by the truncated form of Eq. (jSj) 

N 

fN{x) = ^ + ajTj{x), (9) 

that enables one to evaluate /at at any point x in the interval [—1, +1] without the need to 
carry out interpolations. A method for obtaining the coefficients aj that does not require to 
evaluate the integrals in Eq. (jH)) is described in Ref. It consists in considering the + 1 
zeros of ^V+i for a = 0, 1, ..N, evaluating the expansion Q at x = for a = 0, 1, ..N and 
thus obtaining a set of + 1 linear equations for the coefficients aj. The matrix involved 
that relates the column vector of the f{C,a) to the vector of the aj has elements formed from 
the values Tj{C,a), with j, a = 0, 1, ..A^.. Details can be found in Ref. P] and in textbooks. 
This is the method used to construct Tables ITTllIVI 

The Chebyshev expansion is particularly suited to obtain the integral In^x') dx' of 
the function /tv without significant loss of accuracy. An expansion of this antiderivative 
function in terms of Chebyshev polynomials 

fx N+l 

Fn{x) = / /jv(x') dx' = J2 bjT,{x). (10) 

has the property that the coefficients bj can be easily obtained in terms of the coefficients 
Qj, by means of a matrix usually denoted as Sl, as is described in textbooks as well as in 
Ref. 0]. The basic reason is that the integral from —1 to x of a particular Tj is given by a 
linear combination of Tj(x) with i < j + 1. For example, J^-^ T2{x') dx' = [T^i^x) — 3Ti(x) — 
2ro(x)]/6, and /_!'^ T^ix') dx' = [T^ix) - 2T2{x) + To{x)]/8. The sum in Eq. ^ should 
rigorously go to the upper limit A^ + 1. However, in numerical calculations the (A^ + l)'th 
term is generally ignored. A similar matrix, called Sr, exists in order to obtain a Chebyshev 
expansion of /at (a;') dx' A numerical verification that the accuracy of the antiderivative 
is of the same order of magnitude as the accuracy of the expansion of the function /^v, again 
for f{x) = exp(x), is shown in the second and third columns of Table ITTl 

The derivatives with respect to x of can also be obtained via Chebyshev expansions, 
but in order to maintain a prescribed accuracy, the truncation value A^ has to be inreased 
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Fat - + e-i 








-0.8 


.25(-ll) 


-.51(-9) 


.12(-8) 


.14(-6) 


-0.6 


.ll(-8) 


.51(-9) 


.10(-8) 


-.81(-7) 


-0.4 


.29(-9) 


-.30(-9) 


-.48(-8) 


.37(-7) 


-0.2 


.27(-9) 


-.23(-9) 


.49(-8) 


.24(-7) 


0.0 


.ll(-8) 


-.55(-9) 


.50(-10) 


-.55(-7) 


0.2 


.37(-9) 


-.24(-9) 


-.52(-8) 


.23(-7) 


0.4 


.20(-9) 


-.32(-9) 


.51(-8) 


.41(-7) 


0.6 


.ll(-8) 


.57(-9) 


-.10(-8) 


-.91(-7) 


0.8 


.21(-10) 


-.58(-9) 


-.15(-8) 


.16(-6) 



TABLE II: Coefficients aj and aj x for the expansion of exp{x) for = 9 
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7 


8 


9 


10 


11 


aj 


.32(-5) 


.20(-6) 


.ll(-7) 


.55(-9) 


.25(-10) 


aj X f 


.16(-3) 


.13(-4) 


.88(-6) 


.55(-7) 


.30(-8) 



TABLE III: Coefficients aj and aj x for the expansion of exp{x) 

accordingly. Call /j^^ = df^/dx, f^^ = d^fN/dx^, etc. One of two methods consists in 
taking the derivatives of the Chebyshev polynomials term by term in Eq. (jHl) 

N 

/JrHx)=5^a,7;(")(x), n = l,2,... (11) 

i=i 

The expressions for rj"^ (x) can be given analytically, and hence Z^"* can be evaluated 
numerically at any point x in [—1, +1]. By taking a derivative of a polynomial of order j, 
the result is a polynomial of order j — 1, whose magnitude is of order j times the original 
polynomial. For example, d'^Tj{x) / dx"^ = [xdTj/dx — j'^Tj] /{I — x"^). That leads one to expect 
that the errors in Table HTl for a derivative of order n are related to the coefficient of the next 
to the last Chebyshev polynomial, (Tat+i) times (A^ + 1)". Table UTTl lists coefficients and 
aj X and by comparing Tables ITTl and UTTl one sees that this expectation is borne out. 
A second method consists in writing a Chebyshev expansion for df/dx 

N-l 

dfN/dx=^ + Y,c,T,{x), (12) 
i=i 

8 
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Fjv - + 


In - e^' 






-0.8 


.29(-14) 


.72(-15) 


.24(-14) 


-.44(-12) 


-0.6 


.36(-15) 


.33(-15) 


.21(-13) 


.12(-12) 


-0.4 


.31(-14) 


.44(-15) 


-.32(-13) 


-.15(-12) 


-0.2 


-.56(-16) 


-.22(-14) 


.29(-13) 


.37(-12) 


0.0 


.33(-14) 


.48(-14) 


-.11(-13) 


-.57(-12) 


0.2 


..11(-15) 


-.40(-14) 


-.21(-13) 


.57(-12) 


0.4 


.24(-14) 


.20(-14) 


.44(-13) 


-.37(-12) 


0.6 


.67(-15) 





-.42(-13) 


.25(-12) 


0.8 


.29(-14) 


.88(-15) 


.19(-13) 


-.69(-12) 



TABLE IV: Same as Table HH for = 13 



and by noting that the expansion coefficients cj are related to the coefficients aj in Eq. 
as follows: c^^i = 2NaN, cm-2 = 2(A^ — l)a„_i, and for j < N — 2, Cj_i = Cj+i + 2jaj. The 
error in df^/dx is approximately equal to the magnitude of cat, that in turn permits one to 
determine the value of from the relation cat = 2{N + l)ajv+i 

In the numerical example given in this section the upper value A^ of the sums in the 
Chebyshev expansions was taken as A^ = 9. However, in the numerical solution of the 
integral equation, as described in the next section, A^ = 15. This leads to accuracies of 
the order of 10~^^, as is discussed in the realistic numerical examples described below. In 
order to demonstrate the rapid gain in accuracy for a small increase in the value of A^, 
we show errors similar to those displayed in Table m for A^ = 13. The accuracy increases 
approximately by four or five orders of magnitude as A^ is increased from 9 to 13. 

Once the coefficients of a Chebyshev expansion Qof a function /Ar(x) are obtained, the 
Fourier components f{x) sm{ax)dx and J^^ f{x) cos{ax)dx of that function can also be 
obtained, as follows. If the coefficients dk of the expansion of the function 

M 

f{x) sin(ax) = ^ dkTk{x) (13) 

k=0 

are known, then the integrals J^^ f{x) sin(ax) dx can be easily obtained by applying the 
matrix Sl described above upon the row vector of the coefficients dk, and remembering that 
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Tfc(l) = 1. In order to obtain the coefficients dk one requires the integral 

in view of Eqs. ((Zj). By using the relation 

2Tk{x)Tj{x) = Tk+j{x) + Tik-j\ix) (15) 

the integrals on the right hand side of Eq. (fTill can be carried out analytically in terms of 
Bessel J functions by using the expression |lj] 

T2n+i{x)^^^M=dx = {-irnJ2n+i{a). (16) 
J-1 \J\-x^ 

For Chebyshev polynomials of even order the above integrals vanish. Similarily, one can 
obtain the coefficients of the Chebyshev expansion of /at (a;) cos(aa;) by making use of 

I' T2„(x)-^^^dx = (-l)"vrJ2„(a) (17) 

In this manner the loss of accuracy in the integrals above that takes place for large values 
of a can be avoided. 

Finally, we remark that the Chebyshev expansions can be used on any interval [a, 6] by 
means of the linear transformation 



2 6+a , , 

X = r — (18) 

b — a b — a 



that maps r G [a, b] into x G [—1,1]. 



III. THE INTEGRAL EQUATION METHOD 

Our method for solving the Lippman-Schwinger equation Q is described below for the 
case of one channel and positive energy. The boundary conditions, and hence the choice of 
the Green's function, is appropriate for a scattering situation. Beyond a large radial distance 
called T the potential other than the centripetal or Coulomb potentials is set to zero. The 
radial interval [0,T] is partitioned into subintervals i, with i = 1,2, ...M. The lower and 
upper boundaries of interval i are and bi, respectively, with 6m = T. In each partition 
the integral operator /Cj is defined 

I r 1 f^' 

fCi = — — cos(A;r) / dr' sin(fcr') V{r') — — sin(fcr) / dr' cos{kr') V{r'), 6j_i <r <bi. 

^ Jbi-i ^ Jr 

(19) 
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This operator is similar to ICt defined in Eq. with the exception that the upper and 
lower limits of the integration are and hi. Two independent local solutions Yi{r) and 
Zi{r) in partition i are obtained by solving the integral equation locally, driven by two 
different functions sin(A;r) and cos(A;r), 



It is important to note that boundary conditions are not needed to make the solutions 
of Eqs. fj20|l unique, unless the operator (1 — /Cj) has zero eigenvalues. This situation is 
of course different from the solutions of the differential equation (P), since the functions 
sin(/cr) and cos(A;r) are eigenvectors of the operator {d? / dr"^ + k"^) corresponding to zero 
eigenvalue. If accidentally the operator (1— /Cj) has a zero eigenvalue in a particular partition, 
then by decreasing the size of the partition the zero eigenvalue should disappear because 
the "size" of /Cj decreases correspondingly. Another advantage of the integral equation 
method over the differential equation method is that the operator /Cj is compact, while the 
operator (rf^/rfr^ + A;^) is not. A compact operator can be approximated to ever increasing 
accuracy by a separable expansion of basis vectors, and hence a numerical representation 
(or discretization) of the operator is numerically stable. 

The values of the functions Y and Z and their derivatives at the boundary points of the 
partition i can be obtained from Eqs. (j^UI) by inserting into Eq. ()19|) for r the value or 
6j, respectively. By defining the dimensionless quantities 



{GZ)i = - / cos{kr)V{r)Zi{r)dr ; {FZ)i = - / sm{kr)V{r)Z,{r)dr (21) 



(1 — ICi)Yi =sin(A;r); < r < bi 



(1 - ICi)Zi =cos{kr); bi^i <r <hi. 



(20) 




one obtains 



F,(6,_i) =sin(A;6,_i)[l-(Gr),] 
y/(6,_i) =A;cos(A;6,_i)[l-(GF),] 



Ziipi-i) = cos(A;6i_i) - sm{kbi^i){GZ)i\ 



k[sm{kbi_i) + cos{kbi^i){GZ)i] 



(22) 
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and 



Y,{bi) = sinikbi) - cos{kbi){FY)i 

= k[cos{kbi) + sm{kbi){FY)i] 
Z,ik) = cosikbi)[l-iFZ)i] 
Z[{bi-i) = -A;sin(A;6,_i)[l - {FZ)i]. (23) 

In the above, a prime indicates a derivative with respect to r. Since the functions Y and Z 
obey the Schrodinger equation the wronskian of these functions, Z) = Y'Z — YZ', 

is independent of the point r within the interval i if is a local potential. Using the Eqs. 
(I22|) and (j23|) one can express the wronskian at r = and r = 6j, respectively, in terms 
of the overlap integrals defined in Eq. (|21|). One obtains 

WiY,Z),^_, = k[l~{GYU 
W{Y,Z\ = k[l-{FZ)l (24) 

which implies in particular that 

[GY], = (FZ),. (25) 

This result also shows that if (GY) becomes close to unity in a particular partition, then 
the functions Y and Z will no longer be significantly linearly independent of each other, 
and the JEM method becomes unreliable in this partition. The remedy is to decrease the 
length of the partition, since the value of (GY) will then also decrease. 

The solution of Eqs. (j^Uj) in each interval i is accomplished by expanding these functions 
in terms of Chebyshev Polynomials, and solving the matrix equations for the corresponding 
coefficients. The procedure is well described in Ref. j^l, and will not be repeated here. 
However, a few remarks are in order: 1. The coefficients of the expansion of the functions 
Yi{r) and Zi{r) in terms of the Chebyshev polynomials are obtained with high spectral 



accuracy by using Che 
Clenshaw quadrature 



jyshev collocation points in each partition, together with the Curtis- 
13[. 2. The Eqs. ()20|) are not the inverse of the Schrodinger Eq., 
otherwise there would be no gain in accuracy in using the integral equation. 3. The inverse 
of the operator (1 — /Cj) always exists if the partition i is made small enough, because then 
the operator /Cj becomes small in comparison to the unit operator 1. 4. The calculation 
of the functions Yi{r) and Zi{r) is not computationally expensive, because the number of 

12 



collocation points in each partition is prescribed to be small (16, usually), and hence the 
matrices involved, although not sparse, are of small size (e.g. 16 x 16). 5. The accuracy of 
the calculation of the functions Yi{r) and Zi{r) can be prescribed ahead of time by examining 
the magnitude of the last three coefficients of the expansions. If they are not smaller than 
the prescribed accuracy, then the size of the partition is reduced by a a factor of two, and 
the accuracy will increase correspondingly. This adjustment of partition sizes can be done 
automatically, as is demonstrated in detail in Ref jl5| . 

Next the calculation of the global function ilj{r) in each partition i is described. Since 
the functions Yi{r) and Zi{r) are linearly independent solutions of the Schrodinger equation 
([Q), and since the latter is a linear equation, the function ip{r) can be expressed as a linear 
combination of these two functions 

^(r) = AiY,{r) + BiZ,{r), < r < k. (26) 

A relationship between the coefficients A and B in one particular partition i and those in 
the other partitions can be obtained by returning to the original Lippman-Schwinger Eq. 
(12)) for the function '?/'(r), with r contained in that particular partition i. By expressing the 
integrals in Eq. as sums over the integrals over all partitions, by inserting for ?/'(r) the 
expression for every partition, and by making use of Eqs. fOUj) . one obtains 

M 

A = l- J2 [iGy)jAj + iGZ),B,], ^ = 1,2,...M (27) 

j = i+i 

and 

B, = -Y, [(F'Y), A, + (FZ), B,] , ^ = 1, 2, ...M. (28) 

j = 1 

The 1 appears in Eq. ()27|1 and not in Eq. (f^H|) because the "driving term" in Eq. (jSJ 
is sm{kr) and not cos(A;r). When i = 1 then the sum in Eq. (PHj) is set to zero, which 
requires that Bi = 0. That requirement is compatible with the condition that ip{0) = 0, 
since Zi(0) ^ and Yi{0) = 0. 

The equations (j27|) and (j28p can be manipulated in several different ways so as to increase 
the sparseness of the matrices that define the solutions Ai and Bi. One way, described in 
Refs. and is to subtract from each other Eqs. ^I7\ for consecutive values of i, and 
similarly for Eqs. (j2HI)- By defining the column vectors 



= M ; ^ = h c = (29) 
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one obtains 

/ I M 



12 







M21 I M23 
M32 I 



M 



V 



34 



Mm,m-i I 



/ a, \ 

OtM-l 

\ aM ) 



c 
c 

c 



(30) 



where I and are two by two unit and zero matrices, respectively, and where 



(GF), - 1 {GZ)i 




i = 2,3,. .M 



(31) 



and 



M, 



, i = 2,3,..M. 



(32) 





(FF),_i (GZ),_i - 1 

Note that Eq. (jHUj) generally connects the A and i?'s of three contiguous partitions. For 
example, M2iai + 0:2 + M23a3 = C- 

Another way of combining Eqs. (j27l and l^H|) is to first write them into a (2 x 1) column 
form involving the vectors ctj, and subsequently subtracting equations with contiguous i- 
values from each other, however leaving the last equation in its original form. The result is 

Q 

/ Ti -n2 



r2 —^3 
r3 



«3 



V 7i 72 73 •• 7m-i I / 



«M-1 



c 
c 
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(33) 



where 



1 

-{FY), 1-{FZ), 

1 - (GY), -{GZ), 
1 



(34) 



(35) 
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and 



% = I ° ° I • (36) 

(FY), (FZ), 

It is noteworthy that the first M — 1 equations in (fH^ , 

Ti ai = Qi+iai+i, i = l,2,..M-l (37) 

are equivalent to matching the wave function ifj at the end of partition i to at the start 
of partition i + 1. This can be seen by imposing the two conditions ipi{bi) = ipi+ilbi) and 
V^i(^i) — where ipi is the wave function in partition i given by Eq. ()26|) and where 

tp'i is the corresponding derivative. Inserting into Eq. ^lE^ the values of Yi and Zi or their 
derivatives at either the beginning or the end of a partition as given by Eqs. (j22|) or (|23|). 
respectively, one obtains the result 

A = A+,[l - {GY),+i] - B,^i{GZ),^i 
B,+i = -A{FY), + Bi[l-{FZ)4. 

These two equations are equivalent to Eq. (p?7j) . 
By successive applications of Eq. (jHTj) 

ai+i= (fij+i)"^ Fi ai 

one can relate the values of aj, z = 2, 3, ..M, to ai and then use the last of the (jH^ equations 

M-l / ^ \ 

'-ii ai+aM = (38) 

in order to find the value of Ai. It can be shown that Eq. (j38p is compatible with the 
requirement that Bi = 0. 

Several comments are in order, 
a) The "big" matrices in Eqs. (jHHI) or (jHUj) are sparse, and can be solved by Gaussian elimina- 
tion. Since the number of fioating point operations (fiops) is of order M, the computational 
complexity of the S-IEM is comparable to that of the solution of the differential equation. 
This sparseness property results from the semi-separable nature of the integration kernel /C, 
as is shown in Refs. , which however applies only in the configuration representation 

of the Green's function. This part of our procedure also differs substantially from that of 
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Ref. 3. 

b) The scattering boundary conditions can be implemented reliably. This is because the 
Greens function incorporates the asymptotic boundary conditions automatically. However, 
in the coupled channel case for angular momentum numbers L > 0, the coupled equations 
have to be solved as many times as there are open channels because our Green's functions 
are composed of sin(A;r) and cos(/cr), rather than of Riccati-Bessel functions. We show [3] 
that the desired linear combination of the solutions can be obtained without appreciable 
loss of accuracy, since the matrix required in the solution for the coefficients has a condition 
number not much larger than unity. This means that our various solutions are linearly inde- 
pendent to a high degree, contrary to what can be the case with the solution of differential 
equations. 

c) The method is very economical in the total number of mesh-points required in the interval 
[0, T] because in each partition or spectral collocation method requires very few mesh points 
(like in the case of Gauss-Legendre integration as compared to Simpson' integration), and 
the required length of each partition can be easily adjusted to optimal size based on the 
magnitude of the coefficients of the expansion of the functions Y and Z into Chebyshev 
polynomials, as described before. 

d) The calculation can be distributed onto parallel processors. This is because the functions 
Y and Z, as well as the overlap integrals (jTIj) . required for Eqs. (jHHj) or (jHUj) . can be cal- 
culated separately for each partition independently of the other ones. This is an important 
point, since if the number of channels increases, the number of the quantities (pT|l increases 
accordingly. 

Property c) is also important because, due to the small number of total mesh-points, the 
accumulation of machine round-off errors is correspondingly small. In addition, as is well 
known, integration is numerically more stable than differentiation as discussed for example 
in sections 4.4 and 5.2 on pages 203 and 263, respectively, in Ref. [l2|, and is also shown 
in Tables |n] and IIVI Hence the accumulation of the inherent round-off error is smaller for 
the numerical solution of an integral equation than for the numerical solution of differential 
equations. The small accumulation of roundoff errors in comparison to a finite difference 
method is clearly illustrated in Fig. 1 of Ref. |6j], which compares the round off errors in 
the solution of Bessel's equation obtained via the lEM with that of the Numerov method. 
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IV. APPLICATIONS 



The various features of the S-IEM method will now be illustrated by means of examples. 
The spectral property that high accuracy is reached very rapidly (in principle faster than 
any inverse power of the number of mesh-point in a given radial interval) is illustrated for 
the case of the scattering of an electron from an Hydrogen atom. This is a suitable example, 
because the identity between the incoming electron and the electron bound in the atom leads 
to an additional integral term in the Schrodinger equation, if the Pauli exclusion principle is 
implemented via the Hartree-Fock formulation. Rigorously including this term is difficult for 
the conventional finite difference methods, and various techniques were developed for that 
purpose Q|, and additional references can be found in By contrast, in the JEM method 
this additional integral term is easily incorporated without substantial loss of accuracy jioj ]. 
because the integral kernel is semi-separable. A comparison between the S-IEM and a 
conventional NIEM method [2^ is shown in Fig Q. The L = singlet phase shift was 
calculated for the incident momentum k = 0.2 (ao)~^ and T = 50 ao, while the target 
electron was kept in the ground state of the Hydrogen atom. The figure shows that, as the 
number m of partitions is increased, and accordingly the number of mesh-points m x 16, the 
number of stable significant figures in the phase shift increases very rapidly for the S-IEM, 
illustrating the spectral nature of that method. By comparison, for a method employing 
finite difference techniques based on an equi-spaced set of mesh-points, the number of stable 
significant figures increases much more slowly [2^ for solving a very similar integral equation 
non-iteratively by means of the NIEM method. Although it gives a good illustration of the 
numerical accuracy, this example is nevertheless not very realistic physically because the 
virtual excitations of the bound electron to the myriad of possible states, both bound and 
in the continuum, is not included. Inclusio n o f these excitations requires "state of the art" 
calculations that are presently in progress j2l| . 

Another example is the scattering of atoms at ultra-low temperature. This information is 
needed for the investigation of photo association j^l of the two atoms into a molecule, and 
also in the formation of Bose- Einstein condensates (BE) 2^ . The lifetime of a BE condensate 
is reduced j^l by the three-body process in which two of the atoms combine to form a 
molecule in the presence of a third atom, that in turn carries away the energy of formation of 
the dimer. The depletion rate is proportional to the fourth power of the scattering length. At 
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FIG. 2: Comparison of the numerical stability of two methods for calculating the singlet phase 
shift for electron- hydrogen scattering, as described in the text. The number of significant figures 
on the y-axis is the number of decimal places for which the result remains the same as the number 
of meshpoints is increased. S-IEM is the the spectral method described in this paper, and NIEM 
is a non-iterative method of solving the same integral equation carried out by Sams and Kouri. 

low energies a stable method of calculation is required because, the lower the incident energy, 
the more the long-range part of the potentials contributes significantly to the phase shift. A 
bench mark calculation was performed using the S-IEM method, involving two channels, one 
closed and one open The numerical stability of the L = scattering phase shift as a 

function of the number of mesh points used was investigated, and was compared with various 
other methods of calculation, and the results are shown in Fig. El In all of these calculations 
the maximum radius is T = 500 atomic units (ao or Bohr), the diagonal potentials are of 
the Lenard Jones form C^/r^ + Cu/r^'^, and the coupling between the two channels is of 
an exponential form [l^. At small distances, due to the large depth of the potentials, the 
wave function oscillates rapidly, and hence it is important to be able to adjust the size of the 
partitions accordingly. Since no analytical exact comparison values exist, the "error" in the 
figure is defined as the absolute value of the difference between the result for a given value 
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FIG. 3: Comparison of errors for various methods of computation of the L = phase shift for cold 
atom colhsion, as a function of the number of mesh points in a fixed radial interval. lEM is the 
method described here, FEM is a finite element method, Gordon and LD (logarithmic derivative) 
are two finite difference methods, as explained in the text. 



of the number of mesh points and the maximum va^ 



ue of N employed in the particular 



3i 



method. The FEM method is a finite element method implemented by B. D. Esry and 
carried out by J. P. Burke , Jr [l5|; the Gordon method j26[ was implemented by F. E. 
Mies llSl . and LD is a logarithmic derivative method implemented by the code MOLSCAT 
2^ , j^^- For the LD curve the roundoff errors apparently overwhelm the truncation errors 
when the number of mesh points is larger than 2 x 10^. The S-IEM again shows a rapid 
improvement of accuracy with the number of mesh-points, and it reaches a somewhat higher 
stability than the FEM. Our bench mark calculation was recently used for comparison 
with a finite difference method in which the potential in each partition is assumed constant 
(similar to what is the case with one form of the Gordon method), and the corrections are 
taken into account iteratively. 

In many quantum mechanical calculations, penetration of the wave function through a 
barrier is involved. Examples in nuclear physics are the alpha particle decay of a nucleus, or 
the fission of a nucleus into two daughter nuclei, or in the scattering of a nucleus by another 
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FIG. 4: Numerical error in the phase shift for scattering from a Morse potential with a barrier, in 
the region of a narrow resonance, as a function of the incident momentum. The error is obtained by 
comparison with the analytic result, the momentum closest to the resonance occurs for k = 1.50716. 

nucleus, and also in many similar situations in atomic physics. A barrier frequently occurs 
when a long range repulsive potential, such as a centripetal potential of the form L(L + l)/r^ 
or that of a repulsive Coulomb potential, is added to an attractive nuclear or atomic potential 
of a shorter range. For the scatterin g or the fusion reaction of a light nucleus with a heavy 
nucleus at low incident energies js^, [s^ the penetration of the corresponding wave function 
through such a barrier can pose substantial calculational challenges j32| . In low temperature 
atom- molecule scattering, similar barrier penetration effects become crucial js^. For this 
reason a test of the accuracy of a calculation for a case involving barrier penetration was 



performed. The potential chosen is an " inverted" form of the Morse potential 



for which 



analytic results exist for the scattering phase shift j35|. It has an attractive negative valley 
near the origin at r = followed by a smooth positive energy barrier, a situation which leads 
to resonances. For resonant energies the wave function in the valley region can become very 
large if the width of the resonance is sufficiently small, and in the barrier region this wave 
function decreases as a function of distance. This decrease of the wave function in the barrier 
region amplifies the numerical errors, since in this region the numerical errors tend to increase 
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exponentially. The accuracy of three methods of calculation for a particular resonance 
which occurs for an incident momentum k in the region 1.5071/m~^ < k < 
are illustrated in Fig. |3 The parameters of the Morse potential are given in Fig. 10 of Ref. 
Issj l , the maximum amplitude of the wave function in the valley region at the resonance near 
k = 1.50716/m, is close to 300 (asymptotically it is equal to 1). The error is defined as the 
difference between the analytical and the numerical results; the momenta k on the x-axis are 
given as the excess over the momentum at the left side of the resonance, k = 1.50710 fm~^. 
The lEM curve is obtained with the method described in this paper, NUM is a sixth order 
Numerov method, also denoted as Milne's method Q], and the LP curve is obtained with the 
Logarithmic Derivative method, implemented by MOLSCAT 37]. The matching radius for 
the two finite difference methods, LD and NUM, was set at 50 fm, and the corresponding 
analytical values were extrapolated from T = oo to T = 50 fm by a Green's function 
iteration procedure described in Ref. |[l5|, and are listed in Table 1 of Ref. 3^. For the 
more precise S-IEM calculation that extrapolation was not accurate enough, and T = 100 
was used instead. One sees from the figure that the accuracy of the S-IEM is several order 
of magnitudes (six) higher than that of the Numerov method. 

V. DISCUSSION AND CONCLUSIONS 

A simple way to distinguish a spectral method from a finite difference method is that, 
in a particular partition, the mesh points in the former are not equi-spaced, while in the 
latter they are. Even though the accuracy of finite difference methods can be substantially 
increased by extrapolating the algorithms to equivalent zero-sized distance between mesh 
points js^, such extrapolation methods may become cumbersome. Our spectral S-IEM 
method is one of a class of well-known methods that divide the spatial domain into par- 
titions (or sectors), and expand the solution on a suitable set of basis functions in each 
partition. One example is the method of Gordon [26], that uses Airy basis functions. The 
potential in each partition is approximated by a linear function, and the Airy functions are 
the corresponding exact solutions of the differential equation. This method was included 
among the comparisons carried out for the atom-atom scattering case, illustrated in Fig. 01 
Gordon's method is widely used for atomic physics calculations, and one of the implemen- 



tations can be found in Refs. 



33] and 



291 ]. This is a "potential following method" that 
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is particularly efficient when the potential varies slowly with distance. Another example is 
the method utilized by Light and Walker in which the potential in each partition is 
approximated by a constant. In this case the Green's function that propagates the solution 
from one end of the partition to the other can be written simply in terms of sine and cosine 
functions. This method lends itself well to propagate the inverse of the logarithmic deriva- 
tive of the solution from one end of a partition to the other end, without calculating the 
solution itself. This is called the R-matrix propagation method, and has been implemented 



by Burke and Noble |4l|. This method, as implemented by the code MOLSCAT j27|, was 
included among the comparisons carried out for the barrier penetration calculation, illus- 
trated in Fig. ^ A "function following" method that expands the Greens function in a 
given partition in terms of Legendre Polynomials, without making approximations on the 
potentials, is given by Baluja et al. 



code FARM 



42| . This method is also implemented in the computer 
4l| . The resulting expansion of the distorted Green's function Q (r, r') is of a 
separable form, i.e., it is given as a sum over products of functions uir) x v(r').A similar 

n n 

form is obtained by using Sturmian basis functions |42l], |4^. However such expansions do 
not converge to high accuracy because the derivative of a Green's function has a disconti- 
nuity at the points r = r'. Our S-IEM method does not suffer from that difficulty because 
the distorted Green's function Q (r, r') is obtained in terms of the exact undistorted Green's 
function Qoir^r') through Eq. (plj) . The numerical solution of Eq. is equivalent to 

expressing the distorted Green's function in terms of the undistorted one, according to 

and since Qq is given exactly in terms of its semi-separable form [near Eq. (j2))] there is no 
loss of accuracy. The functions Y{r) and Z{r) are two independent solutions of both the 
Schrodinger equation and the Lippman-Schwinger equation in a particular partition, and 
they represent the two basis functions in terms of which the global solution is obtained in 
each partition. The equation ()37|) . based on algebraic matrix Eq. ()33|) . that relates the 
two expansion coefficients in one partition to the coefficients of one adjoining partition is 
equivalent to the propagation of the logarithmic derivative from one partition to the next. 
However, the method represented by Eq. ()3()|1 relates the coefficients in one partition to 
those in two other partitions appears not to be as closely related to the propagation of the 
logarithmic derivative, hence a comparison of the two methods for particular cases would 
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be very desirable. The method involving two adjoining partitions can be shown to be 
very similar to the multiple shooting method for solving two-point boundary vaue problems 



45[. How the computational complexity scales with the number of coupled channels, in 



comparison with that of other methods, has also yet to be investigated. 

In summary, a recently developed method for solving the Lippman-Schwinger integral 
equation is described and is applied to the solution of several physical problems. Since the 
new S-IEM is considerably more stable than finite difference methods, it is concluded that 
the S-IEM may become the method of choice for particular applications, such as atomic 
physics calculations that involving large distances, require high accuracy, and need to be 
carried out in configuration space. 
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